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1. Introduction 

In evolutionary biology, populations are subject to a number of forces that shape their 
genetic composition [1]. Amongst these, mutations, selection and drift play a central 
role. Drift becomes dominant for small populations, while for large populations one 
reaches a steady state where mutations balance effects of selection. The landscape 
paradigm [2, 3] provides a relation between genotype/phenotype and fitness, allowing 
for quantitative studies of evolving populations, while at the same time giving a 
qualitative picture. This has been particularly developed in the context of quasispecies 
theory [4, 5, 6] and of evolution on networks [7, 8]; in a different framework these 
problems have been analyzed within evolutionary game theory [9, 10, 11, 12] although, 
in that case, usually there is no drift in the usual sense and the number of genotypes is 
low. 

Let us consider a population of M individuals evolving in a static fitness landscape. 
We can define its steady-state distribution from the average number of individuals 
with a given genotype, averaged over a long stretch of time. One can also consider 
how the population is spread out in genotype space. If all genotypes have the same 
fitness (flat landscape), the steady state distribution is independent of M while the 
spread of the population depends mainly on the product M/j, where \x is the mutation 
rate [13]. In the limit of an infinite population (// fixed), the generation-to-generation 
fluctuations vanish and the instantaneous population distribution coincides with the 
steady-state one. This is the quasispecies limit, where the steady-state distribution is 
given by the leading eigenvector of the evolution operator [14]. When the population 
is finite, no general analytic solution for the steady-state distribution is known. If 
the population size is much greater than the number of genotypes, then the so-called 
diffusion approximation [15] can be used. However, in most realistic cases, the number 
of possible genotypes is much greater than the population size. We are interested in 
understanding how mutation, selection and drift affect different properties of the steady 
state in such systems. 

To unravel the different effects, we consider different evolutionary dynamics, in 
which some of these processes may be present or not. Most of our study is conducted in 
the framework of RNA neutral networks [16, 8], the archetypes of genotype to phenotype 
mappings, but we also consider toy fitness landscapes. After specifying our systems and 
population dynamics in sections 2 to 4, we examine the dependence on the mutation 
rate [i and on M of the different processes driving the dynamics. We begin with the 
case of neutral evolution in section 5, turning on and off the drift. In section 6 we allow 
for selection in the usual way that leads to significant drift. As a general rule, drift in 
these situations leads to M/i scaling, as shown previously in the absence of selection [13]. 
Then in section 7 we introduce a particular dynamics with selection but low drift. There 
the Mfi scaling is replaced by a smooth large M limit even when \i — > 0; furthermore 
the genotypic diversity becomes proportional to M rather than being nearly constant. 
The corrections to these scalings generically go as 0(1/M), with large effects when there 
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are barriers in the landscape, as we exhibit in section 8. Section 9 is devoted to some 
final considerations. 

2. RNA neutral networks as fitness landscapes 

In studies of genotype to phenotype mappings, one often focuses on biological molecules 
because the corresponding mapping is relatively well defined. The genotype is simply the 
sequence of the bio-molecule, while the phenotype is its shape, as specified for instance 
by the minimum free-energy structure it folds into. Using either protein or secondary 
RNA structures, it has been found [8, 16, 17, 18] that neutral genotypes (genotypes 
that have a given phenotype) which are connected via single mutational steps form 
extended networks that permeate large regions of genotype space. These are known as 
"neutral networks" . Via the neutral network, a population can move in genotype space 
without crossing unfavorable low-fitness regions, in contrast to what happens in many 
rugged fitness landscapes [2, 19, 20, 21]. However, because of the huge dimensionality 
of our genotype space, large neutral (or nearly neutral) networks can be argued to be 
inevitable [3]. 

Here we shall work with an RNA neutral network, i.e., all RNA sequences which 
fold into a given target RNA secondary structure. The genotype of an RNA molecule is 
given by its base sequence: there are four bases, A,C,G, and U, and thus 4 L genotypes for 
molecules of L bases. The molecule's phenotype is given by its secondary structure, i.e., 
by which bases are paired with which as occur in its folded form. To every genotype one 
associates just one phenotype (the secondary structure of minimum free energy) while in 
general there will be many genotypes compatible with a given phenotype. This many-to- 
one genotype to phenotype mapping has been widely studied [8, 14, 16, 18, 22, 23, 24]. 
Standard computational tools are available on the web to fold given sequences; see for 
instance the fold subroutine from the Vienna package [25] which we used for all of this 
work's computations. Two sequences are nearest neighbors (connected on the neutral 
network) if and only if they differ by a single nucleotide substitution. In general, RNA 
neutral networks are heterogenous graphs, so that for instance the local connectivity 
varies quite a lot from site to site. 

The secondary structure (phenotype) chosen in this study is the one used by 
van Nimwegen et al. [14, 21]: it has 18 nucleotides with six base pairs and is depicted in 
figure 1. By single- nucleotide substitutions, purine- pyrimidine base pairs (G-C, G-U, 
A-U) can mutate into each other, but not into pyrimidine-purine (C-G, U-G, U-A) 
base pairs. Hence we considered only the purine-pyrimidine base pairs. Given the base 
pairing rules for this system, the number of a priori "compatible" sequences for such a 
structure is 4 6 x3 6 = 2985884. At 30°C, 37963 of these fold into the target structure; this 
number depends a bit on the choice of temperature since the fold algorithm computes 
free energies. We find these genotypes to be organized into three neutral networks 
(connected components), of sizes 489, 5784 and 31484 respectively. This will allow us 
to investigate the effect of neutral network size on our observables. 
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Figure 1. The target RNA secondary structure (cf. van Nimwegen et al. [14]). 



To define the fitness landscape, we consider that individuals with the "correct" 
phenotype (residing on the neutral network) are viable, i.e., have maximal fitness, while 
all other phenotypes are non- viable, i.e., have minimal fitness. We take these extreme 
values to be 1 and 0, corresponding to the strong selection limit; then any mutation 
that takes one off the neutral network is lethal. 

3. Population spread, Hamming distances and neutrality 

We now ask how drift and selection affect observables associated with the steady state 
population. A first observable quantifies how much the population is "spread out" in 
genotype space, namely how individuals at a given generation differ from one another 
in genotype space. One defines the "Hamming distance" h between two genotypes as 
the number of positions where the two associated sequences have different nucleotides. 
Following Derrida and Peliti [13], we shall study the distribution P(h) of distances when 
two genotypes are taken at random in the population, averaged over generations. If P(h) 
is broad, then the population is spread out in genotype space. 

A second observable is "neutrality" [26]. Let g be a genotype (a sequence of 
L bases) belonging to the neutral network; examine its 3L possible single-nucleotide 
substitutions and let d be the number of these mutants that belong to the neutral 
network. The "neutrality" of go is then d, the coordination (degree) of go on the neutral 
network. A related notion is the mutational robustness of go- It is defined as the 
survival probability of its mutant offspring. In the context of neutral networks with 
fitness values and 1, we see that in fact 



where d is the neutrality of go. These definitions can be straightforwardly extended 
to the neutrality or robustness of any collection of genotypes. Thus one defines the 
"network neutrality" of a neutral network as the mean of d when considering all of its 
nodes. Similarly, when one has a population of genotypes, the "population neutrality" 
is simply the average of d over that population, each individual being counted once. 
The population neutrality depends on both the neutral network properties and on the 
evolutionary dynamics [14]. Furthermore, we immediately see that the neutrality of a 
population is 3L times the mean robustness of its individuals. 




(1) 
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4. The model 

Many processes affect the genetic makeup of natural populations. In this work, we focus 
on the effects of mutation, selection and drift. We wish to turn on and off selection or 
even drift while considering the effects of the population size M or of the mutation rate 
fj,. Thus the details of the evolutionary process we consider are tailored to emphasize 
one or another of these aspects at a time. 

We consider a population evolving with nonoverlapping generations; the population 
size is kept at a fixed value in the standard way. Other choices could have been made, 
but as in most studies, the detailed procedures used in the evolutionary dynamics is not 
expected to be very important. 

Given the M individuals at the current generation, we must produce M viable 
offspring to form the next generation. Each offspring is produced from a parent and 
given a chance to mutate: with probability 1 — /i, no mutation is applied, and with 
probability \i one base at random is changed. Then selection is applied: the child is 
kept if and only if it is viable. More generally, on an arbitrary fitness landscape, we let 
it survive stochastically according to its fitness. Of course if there is no selection, the 
offspring is always kept. The process of producing offspring is repeated until the new 
generation has size M. Drift comes in via the way the parents are chosen to produce 
candidate offspring. In the standard method, the parents are chosen randomly with 
replacement: clearly this allows for drift as by bad luck some parents will not produce 
any offspring. In the presence of selection, drift cannot be turned off completely but it 
can be significantly lowered. 

Indeed, let us consider the following process. First, each individual of the population 
produces one offspring which mutates with probability ji: if a mutation is lethal, 
the corresponding offspring is killed. In this step there is no replacement and the 
resulting offspring population size will generally be smaller than M. Second, one chooses 
individuals randomly from this offspring population and replicates them. This is done 
until the population size reaches M again. Note that when /j, is small, the new generation 
will be nearly identical to the previous one, even for small M, so there is very little drift. 
Because of selection, a small amount of drift does occur, but its intensity is proportional 
to /i. 

In all our runs, we initialize arbitrarily the population and let it evolve for a large 
number of generations until initial conditions are forgotten: this is the steady state 
limit. All the data presented in this work are time averages taken from this regime. We 
are now ready to see how M, fi, drift and selection affect the spread, neutrality and the 
genotypic diversity of the steady state population. 

5. Dynamics without selection 

In this section, we consider a population evolving without being subject to selection. 
We investigate the effects of allowing or not drift, first on regular networks and then on 
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heterogeneous ones. 

5.1. Homogeneous networks 

Consider the space of sequences of length L; if we take these to all be viable, then we 
get a homogeneous network in which all 1-mutant neighbors of each genotype belong 
to the network. Derrida and Peliti [13] studied the evolution of a population on such 
a flat fitness landscape. Using the fact that there is no selection, it is easy to show 
that the steady state distribution is uniform. Thus the population neutrality is trivial, 
being given by the degree of the network, i.e., 3L, for all population sizes. In contrast 
to this simple result, the distribution of Hamming distances between genotypes in the 
population is generally non-trivial. Depending on the nature of the dynamics, we have 
the following behaviors: 

(a) Drift off — The offspring are produced from parents without replacement: since each 
individual has exactly one offspring, each lineage acts like an independent random 
walk. Thus the Hamming distances between the genotypes in the population 
are completely random: the mean of h lies at 3L/4 (at each position along the 
sequence, one has a 3/4 chance of having different bases when comparing two 
random sequences) and its variance is equal to 3L/16. 

(b) Drift on — Here the offspring are produced from parents with replacement: the 
number of offspring of an individual is variable, leading to tree genealogies. This 
situation incorporating drift was studied by Derrida and Peliti [13] and leads to 
a non-trivial P{h) which depends on M and fi. At any given generation, the 
individuals have mutual distances that reflect the fact that they descend from a 
common ancestor, giving rise to a clustering of the population that fluctuates from 
one generation to the next. For our purposes here, we focus on the result [13] that 
the relevant parameter when M is large is Mfi: in particular, P{h) depends only 
on the product M/x at large M, a property that we call M\i scaling. 

5.2. RNA neutral networks 

We now consider a population evolving on an RNA neutral network, defined as the 
subspace of sequences which fold (at 30° C) into the target secondary structure shown 
in figure 1. RNA neutral networks are generally heterogenous. Evolutionary dynamics 
without selection can be implemented by simply "forbidding" attempts to apply lethal 
mutations. There are two natural ways to do this, referred to as blind and myopic ant 
dynamics [27]. In myopic ant dynamics, also called adaptive random walks [19, 20, 28], 
an offspring that mutates is forced to choose a single point mutation that is non-lethal 
(all non-lethal choices are equiprobable) . In blind ant dynamics (also called gradient 
random walks), a point mutation is chosen at random (lethal or not): if it is non-lethal, 
it is accepted, while if it is lethal, it is refused and the offspring is taken to be non- 
mutant. Both the blind and myopic dynamics can be implemented with or without drift, 
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Figure 2. Distribution of Hamming distances h for different values of \i for the myopic 
ant dynamics in the presence of drift. The population size is M — 20. 

according to the method sketched in section 4. Although these dynamical processes 
may appear to be somewhat artificial, they do provide solvable cases. Furthermore, 
van Nimwegen et al. [14] have shown that in the limit of small M/i, the standard 
evolutionary dynamics converges to the blind ant dynamics. 

Consider first the case without drift, in which the sampling of the parents takes 
place without replacement. Then each lineage performs an independent random walk 
on the whole neutral network. As shown by van Nimwegen et al. [14], the steady-state 
distribution for blind ants is uniform on the neutral network, while for myopic ants the 
probability of being at a node of the neutral network that has degree d is a constant times 
d. (Note that if d is the same for all nodes as in regular networks, we obtain the uniform 
distribution as expected.) Given the steady-state distribution, the histogram P(h) of 
Hamming distances is determined from the fact that the lineages are independent. At 
large L, one expects it to become peaked, neutral networks being widely spread out in 
genotype space. Furthermore, the steady-state distribution and the P{h) are M and \i 
independent. 

Let us now allow for drift. The sampling of the parents takes place with replacement. 
Interestingly, the steady-state distribution of the population is not affected by the 
drift: this is due to the fact that the heterogeneity of the neutral network does not 
affect the chances of appearance of an offspring. However, the lineages are no longer 
independent since the population typically shares a recent common ancestor. As 
a consequence, the Hamming-distance distribution P(h) is not determined from the 
steady-state distribution: it is non-trivial and depends on M/i at large M. Since this 
result holds in a far more general context which includes selection, we postpone its proof 
to the next section. Corrections to the M/i scaling are 0(1/M), with typically a rather 
small prefactor. In the M/i — > limit, just as in the general population dynamics [14], 
one recovers the blind ant dynamics. 

For illustration, we show in Figures 2 and 3 the distribution of Hamming distances 
h for a population with myopic ant moves in the presence of drift on a neutral network. 
We see that in spite of the heterogeneity of the network, the M/i scaling holds just as 
for a homogeneous network. 
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Figure 3. Distribution of Hamming distances h for different values of M/i for the 
myopic ant dynamics in the presence of drift. For each M\x we studied M — 20 and 
M = 40. The data with same M/x superimpose perfectly, exhibiting the M[i scaling. 



6. Dynamics with selection and drift 

6.1. The infinite-population (quasi- species) limit 

In the infinite-population limit, M — > oo (with /x fixed), drift is absent and only the 
effects of mutation and selection show up; this is called the quasi-species regime [4, 5, 6]. 
As shown by van Nimwegen et al. [14], the steady-state distribution is given in this limit 
by the dominant eigenvector of a linear operator defined by the adjacency matrix of 
the network. This eigenvector does not depend on /x, but its eigenvalue A does. If each 
individual in the population produced only one offspring, and the unviable ones were 
eliminated without replacement, then the population size would decay by a factor Ao at 
each generation. Relating this decay to the mutational robustness of the steady-state 
population we immediately obtain 



which yields the population's average neutrality (d)^ via equation (1): both and 
(d)oo are \i independent. 

These results can be compared to the population neutrality in the case of 
blind/myopic ant moves, in which there is no explicit fitness-based removal of 
individuals. In the case of blind ants, the probability of residing on any node of the 
network is uniform. The average neutrality seen by such a walker is just equal to the 
average network neutrality. For a myopic ant, the probability of choosing any node on 
the network is proportional to the degree of the node. Neutrality is slightly higher in 
this case and is given by the ratio of second and first moments of the node degrees. 
A standard variational principle [29] shows that the population neutrality is always at 
least as large as the network neutrality, defined as the average degree of the neutral 
network. Of course on a homogenous graph, population neutrality, blind ant neutrality 
and myopic ant neutrality are all equal to the network neutrality. We refer the reader 
to the work of van Nimwegen et al. [14] for a thorough discussion. 




(2) 
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6.2. Case of an infinite population at fixed M/j, 

In biological situations, the mutation rate is very small, making the large-M limit of 
little use: indeed to get to the infinite M limit just discussed, the quantity Mji must 
be large. When M/j, is finite, drift occurs, and it is of interest to consider the large M 
limit at fixed M/j,. In the absence of selection, the M/i scaling law has been derived by 
Derrida and Peliti [13]. In the presence of selection, it has been exhibited from numerical 
simulations by van Nimwegen et al. [14]. It turns out that this scaling follows from the 
dynamical equations, whether or not there is selection: these equations are invariant 
under a simultaneous rescaling of time, M and /i, as long as M/i is fixed, as we now 
show. 

Let Ni(g) denote the number of individuals in the population residing at the neutral 
network node i at generation g. To go to the next generation, let us choose M individuals 
at random with replacement, and let them mutate with probability /i. One then has 

Nt{g + 1) = E Pj(M Pj ) + G t (Mq t , M Qi ), (3) 

(ij) 

where denotes the sum over the nodes j which are nearest neighbors of node %. Also, 
Pj = fiNi/ (3LM) is the probability to choose an individual of genotype j and to have it 
mutate to node i\ Pj(Mpj) is a Poisson random variable of mean Mp 3 \ qi = (1 — //) iVj/M 
is the probability to choose an individual of genotype % and to leave it without mutation; 
Gi(Mqi, Mqi) is the sum of M 0-1 random variables of mean q t , and thus at large M is 
a Gaussian whose mean and variance are both given by Mqi. 

Since we are interested in the large M limit with M/i fixed, let us define Xi = Ni/M 
to be the fraction of individuals residing on node i. If we average equation (3), we recover 
the deterministic (quasi-species) evolution equations for the x^s. But fluctuations do not 
go away at large M if M/i is fixed, instead the intensity of drift goes to a limit. Extracting 
the mean from the Gaussian of equation (3), the stochastic evolution equations for the 
Xi take the form 

Xi (g + 1) = (1 - nMg) + ]T P 3 (M Pj )/M + G t (0, Mq t )/M. (4) 

(ij) 

Summing this expression over M steps we obtain in the limit of large M (and thus 
^-0): 

Axi = Xi(g + M)- Xi (g) = -Mfix t + ^(M^Xj/SL + G t (0, x t ), (5) 

(ij) 

where we have used the fact that the sum of the M Poisson variables contributes via its 
mean but its variance (equal to pf) becomes negligible. Clearly, the steady-state behavior 
of these stochastic equations depends only on the scaling parameter M/i. As expected, 
by averaging these equations, we recover the quasi-species evolution dynamics. Strictly 
speaking, we have shown that there is a limit when evolving M individuals, neglecting 
the decrease in the population size. But restoring the population to a fixed size M 
involves no mutations and so falls into the standard case of drift for a single genotype; 
that stochastic process also reaches a limit at large M when M\x is fixed, so we can 



Population size effects in evolutionary dynamics 



10 




h 

Figure 4. Distribution of Hamming distances at M\x = 0.2 when M = 10, 20, 40 and 
100. The data superpose perfectly. 
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Figure 5. Distribution of Hamming distances at M\x = 2 when M = 10, 20,40 and 
100. The data collapse is excellent when M > 20. 

conclude that the full process (which maintains the population size at M) also depends 
only on the scaling variable M/i as M — > oo. Note that within this scaling framework, 
we recover the \i fixed, M — > oo case of equation (2) by taking Mji — > oo: we shall 
see that the corrections to this limit are linear in l/M/i. One can also consider the 
limit Mfj, — > 0: there, the population structure typically collapses to just one genotype 
at a time, and as shown by van Nimwegen et al. [14], the effective dynamics reduces 
to a random walk on the neutral network, so the population neutrality is given by the 
network neutrality. 

6. 3. Hamming distances in a finite population 

We studied the distribution of Hamming distances between individuals in the steady 
state population evolving on our three RNA neutral networks. Here we report our 
results only for the largest network (of size 31484), as qualitatively similar results were 
obtained with the two other sizes. 

If we fix Mfj,, we obtain M-independent results when M is large, in agreement with 
the scaling law derived in the previous section. However, the value of M at which this 
scaling arises depends on Mfx. As shown in figures 4, 5, 6 and 7, the scaling sets in for 
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Figure 6. Distribution of Hamming distances at Mji = 5 when M = 20,40 and 100. 
Here the data collapse requires M > 40. 
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Figure 7. Distribution of Hamming distances at Mfj, = 10 when M = 20, 40, 100, 200. 
The Mil scaling is good only for M > 100. 



larger and larger values of M as the value of Mfj grows. For instance, when Mfj = 10, 
one needs M > 100 to really see the Mfi scaling convincingly (cf. figure 7). Moreover 
corrections to this scaling go as 0(1/M), i.e., they are of the same form as we found 
in section 5 in the absence of selection. This is a generic property that will be further 
studied in section 8. 

Within the Mfj, scaling, we see the population spread increases monotonically with 
Mfi. In particular, as Mfj, — > 0, the spread goes to 0, while as Mfj — > oo, the population 
spreads across the whole neutral network. 

6.4- Neutrality in a finite population 

We first examine the population neutrality (d) M . We are interested in seeing how large 
M should be for the Mfj scaling to set in, considering in particular the dependence on 
the neutral network size. 

6.4-1- Small neutral network Figures 8(a) and (b) show the average neutrality (d) M as 
a function of M and Mfj respectively, for fj = 0.01, 0.1 and 0.25. 

The average population neutrality depends both on M and Mfj. If fj is fixed and we 
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Figure 8. Population neutrality (<i) M as a function of 1/M and 1/M/i for = 0.01, 0.1 
and 0.25 in our small size network. One can see the M/j, scaling with small corrections 
due to finite M; note also the linear behavior at the origin. 
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Figure 9. 1/M law for the corrections to Mfj, scaling. 



take M — > oo, we recover the quasi-species limit which is //- independent. Similarly, for 
fixed M/i, we obtain the M/i scaling regime by taking large M, just as in van Nimwegen 
et al. [14]. Furthermore, we see that the corrections to the large-M/i limit of population 
neutrality are linear in l/(M/i) (cf. the linear behavior at the origin in Figure 8(b)). 
Finally, the approach to the large M limit at fixed M/i has measurable 1/M corrections. 
(We also find that the value of M/i affects the time taken to approach the steady state.) 
At fixed M/i, the dependence on M of the population neutrality or of the distribution 
of d is rather mild, though more marked than in the myopic or blind ant dynamics. For 
several values of M/i, we show these 1/M corrections in figure 9. From all these data, 
we conclude that the population neutrality has the form 

<d> M = /Mx)(l + ^ + •••), (6) 

where /(M/i) = (d)^ is the M — > oo limit of (d) M at given M/i and 

/(M/i) = /(oc)(l + ] |- + •••), (7) 

describes how the large M/i limit converges to the quasi-species case. Finally, as shown 
by van Nimwegen et al. [14], /(M/i) tends to the network neutrality as M/i — > 0. 
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6.4-2. Medium and large networks Similar results are found for our medium and large 
networks. Figure 10 shows the population neutrality as a function of 1/M/x, exhibiting 
good Mji scaling. However the corrections to this scaling are larger than those found in 
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Figure 10. Average neutrality as a function of 1/M/x for different values of fj, on our 
medium and large size networks. 

the small network; in fact, to see the asymptotic 1/M law for the corrections to the M/j, 
scaling, one has to go to quite large values of M. Finally, as before, the quasi-species 
limit is reached with corrections 0(l/(M(i)). 

Other results that seems to apply generally is that for fixed //, the population 
neutrality (and thus the population robustness) increases with M, but it decreases at 
fixed M/i. Mathematically, this means that A > in equation (6) while B < in 
equation (7). 

7. Dynamics with selection but low drift 

7.1. Framework 

Random drift reduces the population spread and thus delays the approach to the large 
M limit at fixed Lowering the drift would thus allow one to reach the large M limit 
more easily. Furthermore, one would have a higher mutational robustness of the steady- 
state population for a given population size; this higher survival probability suggests that 
biological mechanisms for reducing drift [30] could be selected for in natural populations. 
In this section we study dynamics in which selective pressures are high but the drift is 
particularly low. 

Our dynamics on a neutral network is defined as follows so that drift effects are 
minimized: 

• For a given population of size M and mutation probability per individual /i, we 
have M\i of the individuals in the population hop to any of their 3L neighbors; the 
others are unchanged. If a mutation brings an individual off the neutral network, 
we kill it, otherwise we keep it. 

• We find the number of individuals that were killed and replace them by randomly 
cloning individuals from the remaining population. 
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Note that in these dynamics we perform a random sampling only for a part of the 
population; in fact, in the absence of selection, there is no drift at all. 

To understand the essential difference from the usual dynamics, we use a branching 
process representation of the evolution and the cloning. An individual is represented by 
a point; at a given time there are M points for M individuals (note that as always M 
represents the number of individuals, not the number of different genotypes). From one 
generation (time) to the next, Mjj, points are replicated leading to branchings while the 
other points just proceed without branching. The time evolution generates branching 
processes in which each individual is represented by a point and is connected to its 
parent at the previous time by an edge. Following this branching process, we obtain 
the descendants of a given individual; one can also consider two individuals, follow 
their edges backwards till one reaches their most recent common ancestor. Because of 
the drift, at large enough times the whole population belongs to just one connected 
component. 

Following Derrida and Peliti [13], let us investigate the genealogy of individuals 
going back in generations. For a subset of k individuals in the population, let Wk(t) be 
the probability that their ancestors t generations ago were all different. To calculate 
this quantity, we first determine the probability x^ that k individuals have k distinct 
immediate ancestors (parents). Define q = /x(l — (R^) M ) where (R^) M is the mean 
robustness of the population. Then assuming q is small and neglecting node to node 
fluctuations in neutrality, we have 

Taking the generation-to-generation processes to be independent, one has Wk{t) = x\, 
which can be approximated to leading order in q by 

w k (t) » exp (- ^^^ tj . (9) 

Hence, unlike the random drift case, the time scale is proportional to M/q. We thus 
define a rescaled dimensionless time variable 

r = qt/M. (10) 

The probability that the two individuals shared a common ancestor at most t generations 
back is given by 1 — W2(t) — 1 — exp(— 2r). Thus the probability density p(r) that 
the most recent common ancestor of the two individuals arose between rM/q and 
(r + dr)M/q generations ago is given by 

p(T) = d(i-«, 2 ( T)) = 2exp( _ 2T) (n) 

Therefore the characteristic time for the most recent common ancestor scales as M//z, 
while it scales as M in the usual dynamics with drift of section 6. 

The distribution of times to the most recent common ancestor can be used to 
obtain the Hamming-distance distribution of the steady-state population. Let (f> v (t) 
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Figure 11. Plot of gcnotypic diversity (c/m and Gm) for dynamics with selection and 
drift as function of InM at M[i = 0.4. 



be the probability that two random walkers in genotype space find themselves at a 
Hamming distance v given that they coincided t generations before. We have 

^ = ¥rwrm^vTT^ -^- 2 ^ < 12 > 

where L is the genome length. Then P v , the probability that the Hamming distance 
between two individuals in the steady state population is equal to v, is given by 



/ drp(r)Mt)- (13) 
Jo 

Using the expression of <p v {t) from equation (12), one sees that P v depends on M but 
not much on in particular, one has a well defined /i — > limit at fixed M. This is to 
be contrasted with the random drift case where P v depends mainly on M\i. 

7.2. Consequence for the genotypic diversity 

A high level of genotypic diversity in a population is usually advantageous for survival. 
Here we study how low drift can greatly enhance this diversity by examining two 
measures of the number of different genotypes, namely the actual number qm (which 
is frequency independent), and the inverse participation ratio Gm of the genotypic 
abundances. Explicitly, for a population of size M, if the number of individuals with 
genotype i is rrii, we define 

G„ = (14) 

where the sum runs over all the different genotypes present in the population. 

We show in figure 11 the two different measures of genotypic diversity for dynamics 
with selection and drift. The bottom curve corresponds to Gm, the top one to qm- 
We find that the "absolute" genotypic diversity qm, (taking into account genotypes 
of arbitrarily low frequencies), grows logarithmically (and thus rather slowly) with M 
at fixed Mfi. On the other hand, the rare genotypes contribute less to the inverse 
participation ration Gm, and we find that this measure of genotypic diversity saturates 
in the large M limit at fixed Mji. 
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Figure 12. Plot of genotypic diversity (<?m and Gm) as a function of population size 
M for dynamics with selection but low drift at M/i = 0.4. 




Figure 13. Distribution of Hamming distances for evolutionary dynamics with low 
drift. The distributions depend strongly on M and only very weakly on /i. 



In Fig. 12 we display the same quantities for our low-drift dynamics. We now see 
that both measures of diversity grow linearly with M at fixed Mfi; thus each genotype in 
the population arises just a few times as M — > oo if M/i is fixed. Clearly the reduction 
of drift in this dynamics allows for a high genotypic diversity. 

7.3. Hamming distances 

We can study Hamming distances as was done in section 6.3 for our low-drift dynamics. 
We find that the distribution of Hamming distances depends strongly on M but not 
much on ji and approaches a /i-independent limit as M grows large (figures 13 and 14). 
This is consistent with the above calculation for the scaling laws. As expected, for fixed 
M, there is a non-trivial limit distribution as /i — > 0. 

7.4- Population neutrality 

In the infinite population limit, the value of population neutrality is independent of 
drift effects, and it is also independent of fx. Consider first our small neutral network 
with 489 nodes and network neutrality 10.4499; we determined the largest eigenvalue 
of the adjacency matrix (via Mathematica [31]), obtaining (d)^ = 11.5107. In a finite 




Figure 14. Distribution of Hamming distances for two different values of population 
size (20 and 40) when M/i=0.4, showing the absence of Mfi scaling when the dynamics 
has low drift. 



population we find that (d) M approaches (d)^, with 1/M corrections that do not depend 
much on the mutation rate. In particular, we find 



(d) M = (d) c 



where A = -0.328 ± 0.002 for fx = 0.1 and A = -0.317 ± 0.005 for li 
corresponding fits are shown in figure 15. 



(15) 
0.25. The 



11.512 
11.51 - ' 
11.508 - 
11.506 - 



\y 11.502 - 



11.498 - 

11.4% - 

11.494 - 

11.492 - 
11.49 







0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.01 



1/M 
(li =0. 



1) 



11.512 






11.51 






1 1 .508 


•: 




1 1 .506 






1 1 .504 












11.502 






11.5 






11.498 






11.496 






1 1 .494 






1 1 .492 






11.49 







0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.00-1 0.0045 O.0C 



1/M 

(li =0.25) 



Figure 15. Population neutrality vs 1/M on the small size network. The dotted line 
is a linear fit as in Eq. (15). 

On our medium-size neutral network (with 5784 nodes), we find that the network 
average neutrality is equal to 10.6888, while (d)^ = 12.592 ± 0.0002. We again 
find a 1/M convergence as in Eq. (15) with A = -0.752 ± 0.035 for li = 0.05 and 
A = —0.662 ± 0.02 for li = 0.25 (figure 16). Just as for the small neutral network, A 
does not depend much on ji. 

Finally, on the large neutral network (with 31484 nodes), we find the network 
neutrality to be 12.116, whereas (d)^ = 15.434. Figure 17 confirms the 1/M convergence 
with A = -0.927 ± 0.045 for li = 0.05 and A = -0.944 ± 0.044 for li = 0.25. 

The overall pattern is thus that the data are well represented by equation (15), with 
an A(fi) that grows with increasing network size and depends slowly on /i. We have 
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Figure 16. Population neutrality vs 1/M on the medium size network. The dotted 
line is a linear fit as in Eq. (15). 
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Figure 17. Population neutrality vs 1/M on the large size network. The dotted line 
is a linear fit as in Eq. (15). 
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Figure 18. Population neutrality (d) M vs M for /i = 0.01,0.05,0.1 and 0.25 on the 
large neutral network, showing the insensitivity to /i. 

also checked directly that (d) M is rather insensitive to the value of fi (cf. figure 18). 
Furthermore, in all cases, there is no Mfi scaling (cf. figure 19), in contrast with what 
happens for the case of standard drift (see section 6). 
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Figure 19. Population neutrality (d) M vs Mfi for fi = 0.01,0.05,0.1 and 0.25 on the 
large neutral network, showing the absence of Mjx scaling. 
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Figure 20. The continuous line shows the distribution of the number of neighbors 
for random nodes on the largest network (with 31484 nodes). The curve closest to it 
shows the steady state distribution for two data sets, with /x = 0.1 and 0.25 when the 
population size is 50: they superpose extremely well. The third curve corresponds to 
M = 1000 and fi = 0.1, 0.25 and 0.5, again with excellent superposition. Superposition 
(independence on n) becomes exact when M = oo. 

7. 5. Distribution of neutrality 

As a last point, let us consider the whole distribution Pm(oI) rather than just the average 
population neutrality (d) M . In figure 20 we display several cases of interest for our largest 
neutral network. The left-most curve is for genotypes chosen randomly and uniformly 
from the neutral network. The next curve on the right is for a population of size M = 50 
undergoing selection with low drift; there are in fact two data sets displayed, one for 
fx = 0.1 and one for fx = 0.25. The last curve is for the same algorithm at M — 1000 for 
three values of fx, namely fx — 0.1, 0.25 and 0.5. 

Several comments are in order. First, as M increases, the overall trend is for Pm(^) 
to shift to larger d; this is in agreement with the general property that mutational 
robustness grows with increasing genotypic diversity. Second, there is hardly any 
dependence of these data on fx, a feature particular to dynamics having reduced drift. 
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8. Reaching the large M limit and barriers in the fitness landscape 

8.1. Motivation 

For large populations in the presence of drift, the relevant parameter appears to be M/i, 
as we have seen in sections 5 and 6. When one wishes to evaluate the large-M limit at 
given /i, it is better to minimize drift as in section 7, especially when \x is small, since 
\i plays little role in this modified dynamics. We also saw that the corrections to the 
M = oo limit go like 0(1/M) on neutral networks. Simple arguments [32] suggest that 
this is a generic property of growth with diffusion phenomena, and we shall now confirm 
this on toy landscapes: the 0(1/ M) corrections are not an artefact of the fitness being 
and 1 as we have assumed so far. Furthermore, we wish to get some insight into the 
size of this correction. For our three neutral networks, we found that the factor A of the 
A/M correction grows with increasing neutral network size; however, barriers (entropic 
or fitness) in the landscape are likely to affect A, as we will now illustrate using a few 
toy landscapes. 

8.2. Evolution in a toy fitness landscape 

To build up our intuition, we will consider a space where genotypes are parametrized by 
a real number, and mutations correspond to small changes of this number. Evolution in 
low-dimensional fitness landscapes has been considered by several authors in the recent 
years and has led to a number of insights (see, e.g., [33, 34, 35, 36]). In this case, it is 
convenient to consider a continuous-time limit because this allows for a Fokker-Planck 
formulation. We can start with all individuals at the same position or place them 
randomly in the landscape. After some time, the population reaches a well defined 
steady state. For the numerical simulation of such evolutionary dynamics, we discretize 
time using a time step At; this then gives the following update rules: 

• At each step (small) mutations arise; this means that the genotype x is changed by 
Ax where Ax is a Gaussian random variable of standard deviation \JlDAt [37]. 

• Given the new positions of the individuals, we allow for replication according to 
fitness given by a function —V(x) in the Fokker-Planck language. V(x) is low (or 
even negative) if the genotype has a high fitness; it is large and positive for an unfit 
genotype. Then for an interval At of time, an individual of genotype x will be killed 
with probability 1 — exp(— V(x)At) if V(x) > 0; if instead V(x) < 0, a clonal birth 
will be produced with probability exp(— V(x)At) — 1. Following standard practice, 
if one wants to keep the population size fixed on average near some target value M, 
one simply shifts additively V(x) so that the expected population size is precisely 
M. 

In the M — > oo limit, the details associated with keeping the population at its target 
value no longer matter and the overall process can be formulated as a rate equation. Up 
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to the rescaling of the population to keep its size fixed, the density of genotypes p{x) 
follows the deterministic equation 

^ = i^-V(.WM). (16) 

This is a linear evolution law that is the continuum analog of the quasi-species dynamics. 
At large times the shape of p(x) converges to the eigenfunction of the linear operator 
on the right hand side whose eigenvalue is largest. One can recognize (16) as being an 
(imaginary time) Schrodinger equation. The problem of the steady-state distribution is 
thus mathematically a simple one that can be solved analytically for particular choices 
of the function V(x). Now we can address the question of how the large M limit is 
reached. 



8.3. Harmonic well 

We first consider the case where V(x) = x 2 , which corresponds in the Schrodinger 
equation framework to diffusion in a harmonic well, a case with no barriers. In the 
M — > oo limit, the probability distribution of the population in this landscape is 
P(x) = exp(-x 2 /2)/v^2T. 

For a finite population of size M, we evaluate numerically the steady state 
distribution Pm(x). We see a clear convergence of this distribution to its large-M 
limit, with 1/M corrections: 

P M {x) = P{x) (l + ^\ + O (M- 2 ) . (17) 



The 1/M nature of the convergence clearly appears in figure 21, where we see that 
the amplitude of these corrections is small. This kind of convergence has been justified 
before [32] in the context of population algorithms for solving linear evolution operators. 
Furthermore, it is possible to show that the correction function K(x) goes to a constant 
at large x. When we obtain such a data collapse, we know M is large enough for one to 
extract P(x). 
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Figure 21. (a) Relative deviation P M (x)/P(x) - 1 vs x for M = 6, 8, 10, 12, 15 and 
20. (b) Data collapse plot: M(P M (x)/P(x) - 1) for M = 6, 8, 10, 12, 15 and 20. 
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Figure 22. Data collapse plot for M 1 M 2 (Pm 1 ( x ) - Pm 2 (x))/(M 1 - M 2 ) vs x, for 
several population sizes Mi and M 2 , in the case of a symmetric double well. We have 
taken (M 1 ,M 2 ) = (6, 10), (10, 20), (6, 20), (20, 50) 

8.4- Symmetric double well 

We now consider the case of a landscape with two degenerate optima separated by 
a barrier (passage of low fitness). For that we take V(x) to be an even polynomial of 
degree 4 in x: V(x) = V(—x). Since P(x) is not analytically known, we examine instead 
the quantity MiM 2 (Pm 1 (^) — Pm 2 (x))/(Mi — M 2 ) for different population sizes Mi and 
M 2 . If the convergence goes as 0(1/M), then, provided Mi and M 2 are sufficiently 
large, the data should collapse onto a limit function. This is indeed what we find; the 
case V(x) = x 2 + O.lx 4 is used for illustration in figure 22. 

8.5. Case of an asymmetric double well 

To go from one fitness peak to another, one has to cross a barrier. The previous double 
well has a symmetric steady state, and even with a small population, this symmetry is 
realized. However when the well is asymmetric, finite population effects will be quite 
larger. Indeed, in the M — > oo limit, even a small non symmetric part (V(x) not even 
in x) will lead to a distribution practically concentrated in one well. (This is called the 
"flea and elephant phenomenon" well known in quantum mechanics: when the barrier 
between the two wells is high, even a tiny difference in V between the two sides leads to 
a big effect, just as when a little itching on an elephant's shoulder leads it to put all of 
its weight onto one side.) When the population is finite, this effect is not so evident, and 
the two wells remain nearly symmetrically populated. One has to go to large population 
sizes M in order to come closer to the M — > oo limit. 

To investigate this effect quantitatively, we consider the asymmetric potential 
V(x) = —x 2 + bx 4 + cx with b = 0.05 and c = 0.002, corresponding to an asymmetry of 
roughly 0.1%. We found the 1/M convergence law, but had to go to M 100 to observe 
the data collapse. This is shown in figure 23; note the large scale of the y axis compared 
to the symmetric double well case: the 1/M corrections are much larger here. When 
the barrier height is increased, one needs even larger M values to see the M = oo limit: 
the proper balance of population on each side of the barrier sets in very slowly in M. 
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Figure 23. Data collapse plot for M 1 M 2 {Pm 1 ( x ) - p M 2 (x))/(M 1 - M 2 ) vs x, for 
several population sizes Mi and M 2 , in the case of an asymmetric well. M — 
200,250,300,400,500. 

In landscapes with more than one dimension, there can also be entropic barriers, that 
is passages that are narrow but not of particularly low fitness; such cases are relevant 
for general neutral networks. 

9. Summary and conclusions 

In general, an evolving population undergoes mutation, selection and random drift. In 
this work we quantified the effect of these processes to untangle the different effects, 
using neutral networks and toy fitness landscapes for illustration. In the case of infinite 
populations, there have been many studies. If Mfj, — > oo, (M being the population size 
and \i the mutation rate), drift is absent and one recovers the quasi-species limit [4, 5, 6], 
with results that developed by van Nimwegen et al. [14] in the context of neutral 
networks. For finite M/i (M — > oo), drift effects are important; there, in the absence of 
selection, Derrida and Peliti [13] derived a number of important results. In this work 
we have considered the case M/j, finite with selection, for both M finite and infinite. We 
derived the M/j, scaling even in the presence of selection. The (finite M) corrections to 
this M/j, scaling are 0(1/M), be there selection or not. When M = oo, we find that 
the quasi-species limit is reached via 1/M/i corrections. In all cases, M\x plays the role 
of an effective population size. These laws are summarized in equations (6) and (7). 
We also found that the amplitude of the correction terms showed a slow increase with 
neutral network size. In practice, the M/j, scaling sets in at relatively small values of M. 
Furthermore, at fixed M/i we showed that the genotypic diversity qm of the population 
increases only logarithmically function of M. 

Finally, we considered a dynamics with low drift in section 7. Drift is effectively 
reduced by a factor \x and thus genotypic diversity always grows linearly with population 
size even if M/j, is fixed. One thereby avoids the M/i scaling law, a useful property if 
one wishes to evaluate the large-M limit at small \x. Nevertheless, reaching this limit 
can be seriously hindered by fitness or entropic barriers in the fitness landscape as we 
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saw in section 8. 
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